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Abstract 

Spatial reaction-diffusion models have been employed to describe many emergent phenomena in biological 
systems. The modelling technique most commonly adopted in the literature implements systems of partial 
differential equations (PDEs), which assumes there are sufhcient densities of particles that a continuum 
approximation is valid. However, due to recent advances in computational power, the simulation, and 
therefore postulation, of computationally intensive individual-based models has become a popular way to 
investigate the effects of noise in reaction-diffusion systems in which regions of low copy numbers exist. 

The specific stochastic models with which we shall be concerned in this manuscript are referred 
to as ‘compartment-based’ or ‘on-lattice’. These models are characterised by a discretisation of the 
computational domain into a grid/lattice of ‘compartments’. Within each compartment particles are 
assumed to be well-mixed and are permitted to react with other particles within their compartment or 
to transfer between neighbouring compartments. 

Stochastic models provide accuracy but at the cost of significant computational resources. For models 
which have regions of both low and high concentrations it is often desirable, for reasons of efficiency, to 
employ coupled multi-scale modelling paradigms. 

In this work we develop two hybrid algorithms in which a PDE in one region of the domain is coupled 
to a compartment-based model in the other. Rather than attempting to balance average fluxes, our 
algorithms answer a more fundamental question: ‘how are individual particles transported between the 
vastly different model descriptions?’ First, we present an algorithm derived by carefully re-defining 
the continuous PDE concentration as a probability distribution. Whilst this hrst algorithm shows very 
strong convergence to analytic solutions of test problems, it can be cumbersome to simulate. Our second 
algorithm is a simplihed and more efficient implementation of the first, it is derived in the continuum 
limit over the PDE region alone. We test our hybrid methods for functionality and accuracy in a variety 
of different scenarios by comparing the averaged simulations to analytic solutions of PDEs for mean 
concentrations. 


1 Introduction 

Diffusion of macromolecules and larger entities is the result of a rapid number of collisions with solvent 
molecules. Often, the detailed information about the solvent molecules is not known and the macro¬ 
molecule (particle) appears to move stochastically. The motion of an individual particle can be simulated 
two different ways (Erban and Chapman, 2009); those in which trajectories are defined on continuous 
domains and those in which trajectories are defined by a random walk on a lattice. Stochastic simulation 
of individual particles is appropriate on time-scales which are much larger than the time between solvent 
molecule collisions and, in such circumstances, the particle trajectories may resemble a Wiener process 
(Brownian motion). Coarse-graining the continuous domain into a lattice of compartments, between 
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which particles undergo a random walk, is a simplification which may be appropriate if the lattice spac¬ 
ing is not too large as to impair good spatial resolution. Both off-lattice (Andrews and Bray, 2004; Erban 
and Chapman, 2009; van Zon and ten Wolde, 2005) and on-lattice (Baker et ah, 2010; Drawert et ah, 2010; 
Engblom et ah, 2009; Lampoudi et ah, 2009; Yates and Baker, 2013; Yates et ah, 2012) individual-based 
stochastic simulations have remained popular paradigms for stochastic diffusion simulations. 

It is also common for mathematical models of diffusive processes to treat diffusion using deterministic 
continuous partial differential equations (PDEs) (Murray, 2003). That is, it is assumed that the copy 
numbers of particles are so large that a continuous distribution of particles can be defined which evolves 
deterministically (Berg, 1993). The effect on the time derivative of the concentration distribution of 
particles due to the random motion of individuals is represented by a Laplacian term in the governing 
equation. Other effects, such as particle flow, reactions and other interactions can be incorporated into this 
mathematical framework by the introduction of additional terms in the governing PDEs. The simplicity 
with which diffusion can be represented in these deterministic models has, in part, lead to its popularity 
in modelling complex physical behaviour. This is notably the case in mathematical modelling of biological 
processes (Adomian, 1995; Field and Noyes, 1974; Schnakenberg, 1979) due to the importance of diffusion 
of either chemicals and/or cells in many biological environments. Some deterministic PDE models have 
the distinct advantage over stochastic models of diffusion that they are analytically tractable so that 
the particle dynamics may, in some circumstances, be predicted without numerical simulation (Zauderer, 
1983). In addition, when numerical solutions are required, there exists a plethora of well-established 
techniques at hand (Thomas, 1995). 

Whilst PDE interpretations of the diffusion mechanisms are convenient and popular they lack the 
ability to appropriately represent systems in which stochasticity is inherent (Erban et ah, 2007). For 
example, fluctuations in the number of particles in a system are especially important for systems in 
which one or more of the constituent species is in low abundance (Gillespie et ah, 2013). Such situations 
appear frequently in biological applications due to the small spatial scales that are often involved and can 
lead to dynamics which are significantly different from a deterministic model which may be postulated 
based purely on macroscopic considerations and assuming large copy numbers (Elowitz et ah, 2002; 
Ghosh et ah, 2007; Gonze et ah, 2002; Swain et ah, 2002; Tian and Burrage, 2004; Vilar et ah, 2002). 
The effects of stochasticity have been investigated in spatially extended systems of the type we consider 
in this paper (Baker et ah, 2010; Gillespie et ah, 2013; Woolley et ah, 2011; Yates et ah, 2012). Due to 
the increase in computer capabilities, stochastic simulation techniques for reaction-diffusion systems have 
found popularity and have been implemented in a number of freely available simulation packages (Ander 
et ah, 2004; Drawert et ah, 2012; Hattne et ah, 2005; Ramsey et ah, 2005; Sanft et ah, 2011). 

Individual-based simulations (off- and on-lattice stochastic simulations) of diffusion are capable of 
capturing stochastic effects which PDEs cannot. However, this is at the cost of computational effort 
which becomes an insurmountable barrier very quickly as copy numbers increase (Gillespie et ah, 2013). 
For extended domains, in which concentrations of particle species vary dramatically, a computationally 
feasible model may require localised regions in which diffusion is modelled using a continuous PDE 
(for regions in which concentration is high) and an individual-based simulation (for regions in which 
concentration is low). For example, in modelling calcium-induced calcium release small numbers of 
ions bind to and facilitate the opening of ion channels thus inducing the release of numbers of orders 
of magnitude more calcium ions. An individual approach is appropriate for the initial binding, but 
infeasible for modelling the concentration of ions after the release (Flegg et ah, 2013). Similar multi-scale 
issues occur in the spatio-temporal modelling of cellular signal transduction (Klann et ah, 2012), in actin 
dynamics in filopodia (Erban et ah, 2013) and in the formation of morphogen gradients in the early 
embryo (Bollenbach et ah, 2005; Kicheva et ah, 2007; Tostevin et ah, 2007; Wolpert, 1969). 

The importance of efficient and accurate hybrid algorithms for coupling models of diffusion between 
different regions of space has resulted in a large body of research in the last decade. Most of the 
models which couple the PDE paradigm with compartment-based stochastic simulations involve the use 
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of significant overlapping regions and/or the averaging of particle flux in order to simulate the necessary 
boundary condition for the PDE (Alexander et ah, 2002; Flekkpy et ah, 2001; Moro, 2004; Wagner and 
Flekkoy, 2004). Ferm et al. (2010) also have a method of coupling a PDF to a compartment-based 
simulation which employs operator splitting and r-leaping for the time evolution of the compartment- 
based model where possible. Recently, hybrid models of stochastic diffusion in which on- and off-lattice 
models are coupled between different regions of space have been proposed (Frban et ah, 2013; Flegg et ah, 
2015, 2012; Klann et ah, 2012; Robinson et ah, 2014). Central to these models is the ability to treat 
individual particles as they migrate from one regime to the other. The individual treatment of particles 
allows for the extension of these coupling algorithms to low copy number regimes. This type of coupling 
technique was also recently applied the coupling of a PDF with particles undergoing off-lattice Brownian 
dynamics (Franz et ah, 2012, 2013). 

In this paper we develop hybrid algorithms which couple a PDE in one region of the domain to a 
compartment-based model in the other. This allows particle concentrations to be resolved on a fine scale 
in regions of the domain where an individual-based description is required, but at a much coarser scale 
in regions where particle numbers are large enough to warrant a continuum description. Rather than 
balancing flux at the interface we use a method which is similar to the ghost-cell method (Flegg et ah, 
2015). A defining characteristic of this method is the individual treatment of particles as they cross 
the interface. In our coupling method a small region of the PDE domain adjacent to the compartment- 
based region is allowed to contribute particles to the compartment-based regime as if it were itself a 
compartment. We will refer to this as a pseudo-compartment. The treatment of this pseudo-compartment 
is complicated. From the perspective of the PDE region it contains a continuous probability distribution 
representing the distribution of particles. However, from the perspective of the compartment regime, 
the pseudo-compartment has compartment-like properties. That is, it contains particles which can jump 
into their neighbouring compartment according to standard compartment-based rules for diffusion. The 
pseudo-compartment’s duality allows for particles to behave correctly as they cross individually between 
the two different regimes. When particles cross over the interface a indicator function with the mass of 
a single particle is added to the probability distribution in the region of the pseudo-compartment. 

A one-dimensional graphical representation of this approach is given in Figure 1. 

We present two separate methods (which we refer to as pseudo-compartment methods (PCM)) for 
conceptualising and implementing this coupling scheme. The first method is necessary for the precise cou¬ 
pling of systems with low copy numbers. The second method is a mathematically justified simplification 
of the first for situations in which large particulate copy numbers are present within the PDE domain. 
We will demonstrate these coupling mechanisms in one dimension only but they are easily generalised to 
higher dimensional problems with (hyper)-planar interfaces. 


2 Algorithms 

2.1 Overview of algorithms for deterministic and compartment-based stochas¬ 
tic diffusion 

In this section we will provide an overview of the two different modelling regimes that will be used in our 
hybrid algorithm, highlighting any intricacies in their definition or implementation. 

Compartment-based modelling 

We will denote the region in which we employ the compartment-based regime Vic- In this region we 
partition the domain into K regular compartments of length h. We choose h to be sufficiently small that 
particles in each compartment can be considered to be well-mixed, reacting only with particles in their 
current compartment. We do not keep detailed information about the specific location of each particle. 
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Figure 1. In the left-hand-side of the domain, flp, the distribution of particles is modelled as a PDE. 
On the right-hand-side, Oc, the particle density is modelled with a compartment-based approach. The 
concentration of particles in each compartment, Ci for i = 1,..., K, of flc is represented by a dark blue 
rectangle (numbers of particles can be found by multiplying by the compartment width, h). The 
probability density in Op is represented by the light-blue curve. The two regimes interact through the 
movement of particles from the left-most compartment in Oc, Ci, to the pseudo-compartment, C-i, 
(whose average particle concentration is represented by the height of the green box on the right of Op) 
and vice versa. The height of the green-rectangle is calculated by integrating the PDE solution over the 
region [0, h) and dividing by h. Thus the ‘number of particles’ in the pseudo-compartment can be 
calculated by multiplying the height of the green box by h. 


rather we simply store the number of particles of each type in each of these compartments. Whilst this 
produces a signihcant saving on computational intensity in comparison to Brownian dynamics models (in 
which the position of each individual is stored and/or updated at each step), we lose information about 
the specific location of each of the particles in a compartment and are forced to assume that each particle 
is distributed uniformly throughout the compartment. 

Considering a general system of N species and M chemical reactions, a propensity function, 
i = 1,..., K and j = 1,..., M + 2, is assigned to each reaction in each compartment (including the 
two “diffusive reactions” by which particles leave compartment i). The propensity functions are defined 
such that aij{t)dt is the probability that the reaction occurs in the compartment during the 
infinitesimally small time interval [t, t + dt]. The propensity function for a specific reaction in a particular 
compartment depends on the number of reactant particles available in that compartment (for all but 
zeroth order reactions), the method by which they combine, the kinetic rate constant and (for all but 
first order reactions) the size of the compartment, h (Erban and Chapman, 2009). Similarly, for each 
species, we can characterise the diffusive jumps of particles between compartment i and its two nearest 
neighbouring compartments as a set of first order interactions with propensities 

= AiDA/h'^, ( 1 ) 

for / = M -|- 1, M -I- 2. Here, Da denotes the macroscopic diffusion coefficient of species A and Ai is the 
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copy number of species A in compartment i. Equation (1) is specific for regular square lattices with no 
sensing of neighbouring compartments. For unstructured meshes, the jumping propensity aij between 
two lattice points can be calculated by considering a finite element discretisation of the appropriate 
PDE over the mesh or, in special cases by considering first passage time problems. For details on the 
derivation of transition rates on unstructured meshes the reader is referred to (Engblom et ah, 2009; 
Yates and Baker, 2013; Yates et ah, 2012). Such a discretisation justifies the notation aij (that is, 
a dependence on the destination for the jump). In other scenarios, even on regular lattices, (where 
the particles being modelled are cells, for example) jump rates can depend explicitly on the number of 
particles in neighbouring compartments in order to mimic “volume- filling” or “density-dependent” effects 
(Baker et ah, 2010). Although we do not explicitly deal with these scenarios in this work, it would be 
straightforward to do so. 

There are two predominant mechanisms for simulating the time evolution of the system. In the 
time-driven approach, a fixed time-step At is designated and each reaction j in compartment i occurs 
with probability aijAt. This method is simple to implement, but less efficient than the alternative 
event-driven algorithm. In order to ensure accuracy in the time-driven algorithm the time-step must be 
chosen to be sufficiently small that single particles do not undergo multiple diffusion or reaction events 
in a single step. This can lead to very small choices of At and consequently many time-steps in which 
no reaction/diffusion events occur. In comparison, an event-driven algorithm simulates the evolution of 
the system by choosing the time-step between reaction/diffusion events from an exponential distribution. 
The most commonly cited example of such a simulation scheme is the Gillespie’s direct method (DM) 
(Gillespie, 1977) in which the time for the next reaction of any sort, t, is chosen from an exponential 
distribution with parameter ag (where oq = X]i=i i® the sum of all the propensity functions). 

In practice r is simulated according to the following formula 

r = —In—, (2) 

ao n 

where ri is uniformly distributed in (0,1). The event that takes place r units of time from the current 
time is chosen from a weighted probability proportional to each of the event propensities aij. 

There are efficient adaptations of Gillespie’s original algorithm available depending on context (Gao 
et ah, 2004; Gibson and Bruck, 2000; Li and Petzold, 2006; McCollum et ah, 2006; Yates and Klingbeil, 
2013). These algorithms are mathematically identical and differ only in the computational bookkeeping 
required to implement them. An efficient algorithm designed specifically for spatially extended systems 
is the next sub-volume method (NSM) from Elf and Ehrenberg (2004), and an adaptation of the next 
reaction method (NRM) from (Gibson and Bruck, 2000) (which is itself based on Gillespie’s original first 
reaction method (FRM) (Gillespie, 1976)). For simplicity and efficiency we choose Gillespie’s DM (event- 
driven) in order to evolve particle numbers in Oc, although we note that our algorithms can be trivially 
modified to incorporate the alternative event-driven simulation methods discussed above or time-driven 
evolution. 

PDE-based modelling 

We denote the region in which we solve the PDE as flp. In this region we solve the relevant one¬ 
dimensional PDE numerically, but to a high degree of accuracy. In particular we solve the PDE 

| = + x.np, (3) 

where R{p) represents a set of, as yet unspecified, reactions. 

A zero-flux boundary condition is implemented on the PDE at the interface between the two regions. 
This is because flux of particles over this interface is not governed by a continuous diffusive flow of 
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particles but rather a discrete jumping processes which mimics it. It is important to ensure that there is 
consistency between the methods by which particles transfer back and forth between regimes since any 
imbalance can create significant numerical artefacts. 

In the examples that follow, we discretise the PDE-domain, Qp, into a regular grid with spacing Ax 
where Ax h. In order that the time-step which we take is not restricted by the requirement for stability 
of the PDE solver we employ the 0-method with 0 < 1/2 which is unconditionally stable. In particular 
we choose 0=1/2 (Crank-Nicolson method) which gives a solution which is second order accurate in 
space. Note that the coupling algorithm that follows is independent of the method which is used to solve 
the PDE, providing that method is sufficiently accurate. 

2.2 Low copy number coupling technique: Algorithm 1 

In both hybrid algorithms presented here, it is important to note that, whilst we implemented the Gillespie 
direct method for the compartment-based simulation (event-driven) and a Crank-Nicolson method for 
solving the PDE part of the simulation (both for accuracy reasons), the choice of methodology does not 
affect the results of the coupling (but may affect the accuracy of the individual simulations). 

It is common to think of the PDE description of diffusion as the evolution of a continuous concen¬ 
tration. However, at low concentrations this quantity is not well defined. Here we consider a situation 
in which there may be insufficient particles to satisfy the continuum limit. According to the definition 
of concentration, p{x,t), given Np identical particles in total over the concentration distribution, the 
probability distribution of finding any particular particle, p'{x,t), is given by p'{x,t) = p{x,t)/Np. The 
distribution p'{x, t) is well-defined for low copy numbers and behaves in a similar way to the concentration 
(albeit normalised by Np). In this situation, the continuous distribution which is solved by the PDE 
cannot be interpreted as a concentration. Instead, we consider p'(x,t) as the probability density to find 
each of the Np particles that are within the PDE domain (we denote the number of compartment parti¬ 
cles Nc such that total number of particles is IV = Np -|- Nc). We specify that in the continuum limit, 
as Np —> oo, the PDE should describe the concentration of particles. As such, we define the continuous 
distribution p'{x,t) to be the probability density to find each of the Np PDE-based particles at position 
X, scaled by the number of particles, such that Np {x,t) da; is the expected number of particles 

in an arbitrary subset, w, of Op. 

Our first algorithm (which is appropriate for situations in which the copy numbers are low) progresses 
asynchronously. The compartment-based regime is updated in an event-driven manner whilst the PDE 
domain is updated in a time-driven manner. We use the model design shown in Figure 1 and describe 
here how to couple the transition of particles between the PDE and compartment-based regimes using 
the pseudo-compartment, C-i. 

Particle numbers in region Oc are updated according to jump propensities given in equation (1). 
We extend the definition of these jump propensities to include transitions to and from the pseudo¬ 
compartment, G-i, that lies inside the PDE domain. When a particle is chosen to jump left from 
compartment Ci in He, into the pseudo-compartment, G-i, its mass is incorporated into the PDE 
solution. The expected number of particles within the pseudo-compartment is increased by one by 
adding 1/h to the PDE distribution in [—h, 0). 

In ric, the propensity for a jump from compartment to neighbouring compartment is given by equa¬ 
tion (1). We modify the propensity to jump from the pseudo-compartment, C_i, into its neighbouring 
compartment, Ci, in Oc, in order to take account of the expected number of particles in C-i. The jump 
propensity from C-i to Ci is therefore given by 

a* = (4) 

where A_i = ^p{x,t)dx is the expected number of particles in the pseudo-compartment. For conve¬ 
nience in what follows we denote the normalised number of particles in C_i as a_i = A_i/Np. 
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Jumps from C-i to Ci occur as compartment-based events in the algorithm. This means that in each 
small time interval (t, t -|- dtp) in which the PDE is updated, none of the PDE-based particles can jump 
into flc- Since p{x,t) is the scaled probability distribution describing the probability of finding these 
particles in the PDE regime, the new distribution must be conditioned on the event that none of these 
particles jumped into the compartment-based regime during the PDE update interval (t, t + dtp). After 
running the PDE solver on the distribution p{x,t) a change of Sp{x,t) must be made to reflect the fact 
that none of the remaining Np particles transferred into Ci via C-i. 

In the small time step dtp each of the expected A_i particles inside the pseudo-compartment should 
be capable of jumping over the interface with a probability of J^jump = Ddtp/h^. By definition 


Sp{x, t) = Np {Pr{x, t|no jump) — Pr{x, t)) 
Pr{no iump\x,t)p{x, t) 

= - Pr(nojump) - 

Pr{no iump\x,t)p{x,t) 

= “i- B - , 

1 -^umpO'— 1 H“ Oycitpj 


(5) 


where Pr indicates the probability for a single particle. Pr{no jumpja;, t) depends on whether the particle 
in question is inside or outside the pseudo-compartment, C-p. 


Therefore, 


Pr{no jumpjx, t) 


1 — P- 

jumpi 


X e C_i 


Sp{x,t) 


p{x^ 1 

1 J^umpa—1 

p(x, t)Tjump(c^—1 1) 

1 T^ump^—1 


X i C-x 

X G C-i 


( 6 ) 


(7) 


As dtp —7- 0, J^jump 0 and we hnd that 


lim 

dtp —^0 


Sp{x,t) 

dtp 


Da-i 


p{x,t), 


x i C-i 


Da-i 


p{x,t) 


D . ^ 


X € C-i 


( 8 ) 


The PDE is therefore solved with an additional source term of size p{x,t)Da-i/h‘^, throughout the 
domain and a sink term inside the pseudo-compartment of size —p{x,t)D/h?. These additional terms 
ensure that the probability distribution inside the PDE domain redistributes correctly as a result of not 
finding pseudo-compartment particles available for jumping in a given time step. 

Changes in the compartment-based regime have an impact on the PDE only when particles cross 
the interface. We have already discussed the change to the relatively straightforward indicator function 
addition of density to the PDE when particles cross from Ci to C-i. However, when particles cross from 
C-i to Cl it would be erroneous to simply subtract 1/ft. from the PDE distribution in C_i (reversing 
the operation for the transfer of particles in the opposing direction). This is because we assume that 
each of the particles in the PDE region has the same distribution function p{x,t)/Np. Therefore, taking 
any one of these Np identical particles across the boundary into Dc will result only in a rescaling of this 
distribution. Recognising that we lose a single particle’s worth of mass, the distribution p{x,t) in Dp 
should now represent the probability distribution of Np — 1 identical particles so we can calculate the 
rescaled distribution as 


Pnew {Xj t) 


Np - 1 
Np 


p{x,t). 



A cartoon schematic illustrating the algorithm is given in Figure 2. The algorithm is also outlined in 
detail in Table 1. It may seem odd that particles appear to be allowed to transferred from the entire 
PDF domain whilst being allowed only to transfer back into the more localised C-i. This is restriction 
is predicated on the assumption that all particles in the PDF regime have the same distribution which is 
proportional to the PDF solution, p(x, t). However, when a particle travels from Ci in the compartment- 
based regime into C-i in the PDF regime, its position is known. This breaks the assumption that all 
particles in the PDF are identical. This issue has been documented previously (Franz et ah, 2012). It 
results in an increase in the variance of the number of molecules near the interface. In (Franz et ah, 2012), 
the variance error can be controlled by adding in an ‘overlap’ region in which particles may be in the 
form of either regime (PDF and off-lattice particle-based). The pseudo-compartment in this manuscript 
plays a similar role to this ‘overlap’ region. In the pseudo-compartment, whilst density is updated using 
the PDF, particles are transferred by the evolution of the PDF (into the remainder of the PDF regime) 
and by lattice-based rules (into the compartment-based regime). This complication is unavoidable for 
algorithms of this nature. The variance is also naturally reduced by increasing the number of particles 
(see Algorithm 2). 

In Algorithm 1, changes to the PDF solution are as a result of diffusion, auxiliary source/sink terms 
for changes in conditional probability, events where a particle moves from the pseudo-compartment, C-i, 
to the hrst compartment, Ci, in the compartment-based regime, flc, and events where a particle moves in 
the opposite direction from Ci to C-i. All of these changes may be represented formally in the following 
PDF (in one dimension) 


dp{x, t) ^ d^p{x,t) f N+ 
dt ^\N- 




Da-i(t) ^C-iix) 


D 


. ( 9 ) 


where Is (a;) = 1 if x € i? and zero otherwise, tj are instants when particles jump from compartment 
to PDF, tj are instants when particles jump from PDF to compartment and Np = Np{t~) and Np = 
Np{t^) represent numbers of particles in Hp immediately before and after times in which the number 
discontinuously jumps at respective moments tj (they can be thought of as left and right limits of Np(t) 
respectively). Note that Np = Np — I. Asa common sense check we can integrate over the flp to derive 
an equation for the change in the total number of particles represented by the PDF regime, Np: 

^ =-^(5(t-tj)+^(j(t-tj). (10) 

•J 3 

As expected, the number of particles changes by integer amounts at the discrete times when particles 
enter or leave the domain. 


2.3 High copy number coupling technique: Algorithm 2 

The algorithm presented in the previous subsection is accurate for both high and low copy numbers. 
However, its complexity leads to computational overheads which may prove costly, especially in situations 
in which copy numbers of species are relatively high and a less sophisticated coupling method would suffice. 
For this reason we derive the following coupling method (which we refer to as Algorithm 2) as a limiting 
case of our original algorithm, valid only for the scenario of high copy numbers. 

We predict the benefits of these hybrid methods to be greatest for applications that have regions in 
which there are low copy numbers (requiring individual-based simulation) and other regions of high copy 
numbers (in which a PDF description will be more efficient). In particular the copy numbers at the 
interface should be high in order to justify the use of the PDF there. Arguably if the copy numbers are 
low at the interface then the interface is positioned incorrectly. 
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Figure 2. A cartoon schematic illustrating the most important features of Algorithm 1 in one 
dimension. In to the right of the interface, particles are represented as individuals in compartments 
(dark blue rectangles) with each individual adding 1/h to the density in that compartment. Particles 
can jump to neighbouring compartments within (dark blue block arrows). In Up, to the left of the 
interface, density is represented by a ‘continuous’ function (light blue curve). Particles in the left-most 
compartment in He (adjacent to the interface) can jump into the PDE regime (blue horizontally striped 
arrow) where they are added by means of an indicator function to the PDE solution in the region 
[—/i, 0) - the pseudo-compartment. The PDE does not include information about the specific location of 
individuals and so when a particle jumps from the pseudo-compartment into ^Ic it must first be realised 
in the pseudo-compartment by subtraction of its contribution to the whole PDE followed by an 
instantaneous jump over the interface. Every instance in which a particle does not jump from the 
pseudo-compartment into flci the conditional probability to find particles in the distribution must 
necessarily be adjusted so that more weight lies outside of the pseudo-compartment. This effect results 
in a net source of distribution outside of the pseudo-compartment in the PDE and a balancing net sink 
of distribution inside the pseudo-compartment such that the total number of particles is conserved in 
the absence of reactions. 


In order to derive our simplified algorithm we allow Np —> 00 . This implies that times between 
transfer events (particles crossing from C-i to Ci and vice versa) are low: 0 < tj+i — 1 and 

0 < tj+i — fj "C 1. Let us define a small time, such that 

1. signihcantly smaller than the total number of particles in the pseudo-compartment. 

Satisfying point 1 ensures that individual jump events may be averaged over (f, t St) and satisfying 
point 2 ensures that St is small on the scale of diffusion allowing for PDE solver time discretisation to 
be implemented at a later stage on the majority of non-jumping particles. Both points 1 and 2 may be 
satisfied for sufficiently large Np. Integrating equation (9) from t to t + St, dividing by St and taking the 
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(i) Initialise the time, t = to and set the final time, tfinal- Specify the PDE-update time-step dtp and 
initialise the next PDE time-step to he tp = t + dtp. 

(ii) Initialise the number of particles in each compartment in Oc, Ai for i = 1,. ..,K, and the distri¬ 
bution of probability density in flp, p{x,to), for x G flp. 

(iii) Calculate the propensity functions for diffusion between the compartment regimes as = AiD/h? 
for i = 1.. .K and j = M + \,M + 2 and for reactions as aij ior i = 1... K and j = 1,... ,M 
using the usual mass action kinetics. 

(iv) Calculate the propensity function for diffusion from the pseudo-compartment, C-i, in Dp into the 
adjacent compartment, Ci, in flc- a* = D ^ p(a;, t)dx//i^ (equation (4)). 

(v) Calculate the sum of the propensity functions, oq = <^i,j + ■ 

(vi) Determine the time for the next ‘compartment-based’ event tc = t + r, where r is given by equation 
(2). 

(vii) If tc < tp then the next compartment-based event occurs: 

(a) Determine which event occurs. Each event occurs with probability proportional to its propensity 
function (see Gillespie (1977)). 

(b) If the event corresponds to aij for i = 1... K and j = M + 1,M + 2 then move a particle from 
compartment i in the direction specified by j. If the particle crosses the interface into pseudo¬ 
compartment C-i then add a particle’s worth of mass to the region C-i i.e. p{x,t -I- t) = 
p{x,t) + f-C-i/h. Here Ic_i is the indicator function which takes the value 1 in C-i and 0 
elsewhere. 

(c) If the event corresponds to propensity function a* then place a particle in Ci. To remove this 
particle from the PDE solution, simply rescale the solution i.e. p{x,t + T) = p{x,t){Np — l)/Np. 
Update Np and Nq to reflect the exchange of particles between regimes. 

(d) Update the current time, t = tc- 
(viii) If tp < tc the the PDE regime is updated: 

(a) Update the PDE solution according to the numerical method described in Section 2.1 using 
p{x, t) as the previous value of the solution. 

(b) Add the auxiliary source and sink terms (see equation (8)) to the numerical solution of the 
PDE. 

(c) Update the current time, t = tp and set the time for the next PDE update step to be tp = 
tp + dtp. 

(ix) If t < tfinah return to step (iii). 

Else end. 


Table 1. Algorithm for the coupling of the compartment-based regime with the PDE-based regime in 
the case of low copy numbers (Algorithm 1). 
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limit as —7- 0 we find 
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j2Ht' -tjW > ■ 


(11) 


For each C-i —>■ Ci transfer event we approximate {Np/Np — 1) = —1/iVp Ri —1/iVp, since Np 
is large. It should be noted that in order to obtain equation (11) from equation (9) one requires St 
be sufficiently small that temporally varying dependant variables are constant on (t,t + St), that is 
p{x,t')dt' « p{x,t)St and similarly for a_i(t) and 1/Np{t). The large number of transfer events 
that occur between t and t + St mean that the fluctuations in the number of transfer events are negligible 
and, as such, they can be predicted deterministically using the propensities given in the previous section 
(equations (1) and (4)): 



^(5(t' - tj)dt' 

J 

j 


DStA_i 

h? 

DStAi 


and 


( 12 ) 

(13) 


Substituting equation (12) into equation (11) and recalling that A_i = a-iNp, we find the second and 
third terms on the right-hand side of equation (11) cancel each other. Consequently, we see that changes 
in the PDF as a result of hybridisation with the discrete subdomain (i.e. not including the Laplacian 
diffusion operator) are confined to within the pseudo-compartment. 


dp j^d‘^p{x,t) 

dt ^ dx"^ 


lim 


Die. Ax) , lc-i{x) 
hSt 




(14) 


Substituting equation (12) into the second term on the right hand side of equation (14) we obtain (correct 
to 0{h?) accuracy; as accurate as the compartment-based simulation) 


^ ^ d^p{x,t) 
dt dx'^ 


lim 
<5t—>-0 


hSt 



tj) - '^Sit' -tj)dt'. 

J 


(15) 


In order to find equation (15) it should be noted that 

\c ('x') r> 

lc.Ax)p{x,t) = lc.i{x) {p{xo,t) +px{xQ,t){x - xo) + ■ ■ ■) = —- {A-i + 0{h?)) , (16) 

where xq is at the centre of C-i and px{xo,t) = 0{h) (due to the no flux boundary condition placed on 
the PDE at the interface between the two subdomains). Evaluating the limit St ^ 0 we obtain 


S{t-tj)j. (17) 

If the number of particles within the pseudo-compartment, A.i, becomes large then it will have a small 
relative variance and hence can effectively be assumed to be deterministic (albeit not necessarily integer 
valued). 


dp{x,t) ^^d'^p{x,t) 


^e.rix) 


dt 


dx^ 
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Equation (17) is the motivating representation for Algorithm 2, which is valid for large Np. Transfer 
events from flp and flc are calculated in the same way as Algorithm 1, however, their effect on the 
PDE domain is different. The source and sink terms in Algorithm 1, required for conditional probability 
changes in the PDE distribution as a result of periods of ‘no jump’, are no longer part of Algorithm 2. In 
Algorithm 2, as particles transfer from C-i to Dc particles are withdrawn from the PDE assuming that 
they were from the pseudo-compartment in the same way that they are added when particles transfer 
from from flc to C-i (as indicated in equation (17)). 

Unsurprisingly, for large Np we obtain a near deterministic hybrid simulation in both subdomains. 
Using equations (12) and (13), equation (15) becomes 


dpjx, t) ^ d‘^p{x,t) 

dt dx"^ 




(^1 


A_i). 


(18) 


The last two terms in equation (18) represents a lattice approximation to a free diffusive flux over the 
interface, a desired outcome of Algorithm 2. 

A cartoon schematic illustrating the simplified algorithm for large copy numbers (Algorithm 2) in one 
dimension is given in Figure 3. 

Table 2 describes the algorithm used to interface the two regimes. Note that when removing a pseudo¬ 
particle from the PDE regime in step (vii) (c) we enforce the condition that p{x, t) > \/h for all x G 
C_i to ensure that when we remove an indicator function corresponding to a particle’s worth of mass 
that we do not make the PDE solution negative. In practise this condition will rarely be enforced since 
the PDE region should be used to describe regions of the domain that have sufficiently high density that 
the individual-based model becomes inefficient. We note, however, that with a fixed interface this will 
not always be possible and that this restriction may lead to slight inaccuracies in the solution. 

We have shown that, theoretically. Algorithm 1 simplifies to Algorithm 2 in the case of large Np. In 
the following section, we will demonstrate this empirically on a number of simple test problems. We will 
also demonstrate the strong agreement between the hybrid simulations of both algorithms and known 
exact analytic solutions. 


3 Results 

In this section we assesses the accuracy of our proposed algorithms in simulating three straightforward test 
problems of increasing complexity. In all the example simulations which follow we choose dtp, Nx = 0.005 
and h = 0.05, so that for every compartment there are 10 PDE mesh points. These parameters will be 
systematically varied later in this section in order to determine the error behaviour of the algorithm. 

3.1 Test Problem 1: Uniform particle distribution 

In Test Problem I we demonstrate that our proposed algorithms are capable of maintaining a steady 
state distribution of particles. We consider the diffusion of non-interacting particles on a domain with 
reflective boundaries which will maintain a steady state distribution. We assume that all N particles are 
initially uniformly distributed across the domain and move with diffusion coefficient D = 1/400, in non- 
dimensionalised units. The PDE that we solve in Dp is the diffusion equation with zero flux boundary 
conditions at either end. 

Figure 4 demonstrates that a uniform steady state is maintained by both algorithms. The disparity 
between the PDE solution and the analytic solution close to the interface in Figure 4 (c) is transient and 
oscillates either side of the analytic solution as time progresses (see Supplementary Material (SM) Movies 
Ml and M2 corresponding to Algorithms 1 and 2, respectively). 
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Figure 3. A cartoon schematic illustrating the most important features of Algorithm 2 in one 
dimension. In flc, to the right of the interface, particles are represented as individuals in compartments 
(dark blue rectangles) with each individual adding 1/h to the density in that compartment. Particles 
can jump to neighbouring compartments within flc (dark blue block arrows). In flp, to the left of the 
interface, density is represented by a ‘continuous’ function (light blue curve). Particles in the left-most 
compartment in flc (adjacent to the interface) can jump into the PDE regime (blue horizontally striped 
arrow) where they are added as an indicator function to the PDE solution in the region [—h,0) - the 
pseudo-compartment (green dashed line). Pseudo-particles (hatched green rectangles with dashed 
borders) corresponding to the mass of the PDE in the pseudo-compartment are allowed to jump right 
(green vertically striped arrow) into the first compartment in Qq- When a particle leaves the 
pseudo-compartment an indicator function of mass corresponding to one particle is removed from the 
PDE in that region. 


3.2 Test Problem 2: Particle redistribution 

As a further example of the propriety of our algorithms we consider a situation in which particles are 
initially uniformly distributed in either tip or in flc with no particles in the complementary region. 
We note that the situation in which all the particles are in Qc and hence the concentration in the 
compartment-based region is much higher than that in the PDE region is not the ideal situation in which 
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(i) Initialise the time, t = to and set the final time, tfinal- Specify the PDE-update time-step dtp and 
initialise the next PDE time-step to he tp = t + dtp. 

(ii) Initialise the number of particles in each compartment in Oc, Ai for i = 1,. ..,K, and the distri¬ 
bution of probability density in flp, p{x,to), for x G flp. 

(iii) Calculate the propensity functions for diffusion between the compartment regimes as = AiD/h? 
for i = 1.. .K and j = M + \,M + 2 and for reactions as aij ior i = 1... K and j = 1,... ,M 
using the usual mass action kinetics. 

(iv) Calculate the propensity function for diffusion from the pseudo-compartment, C-i, in Dp into the 
adjacent compartment in flc- a* = D ^ p(a;, t)dx//i^ (equation (4)). 

(v) Calculate the sum of the propensity functions, oq = <^i,j + ■ 

(vi) Determine the time for the next ‘compartment-based’ event tc = t + r, where r is given by equation 
(2). 

(vii) If tc < tp then the next compartment-based event occurs: 

(a) determine which event occurs according to the method described in the text (see Gillespie 
(1977)). 

(b) If the event corresponds to aij for i = 1... K and j = M + 1,M + 2 then move a particle from 
compartment i in the direction specified by j. If the particle crosses the interface into pseudo¬ 
compartment C-i then add a particle’s worth of mass to the region C-i i.e. p{x,t -I- t) = 
p{x,t) + f-C-i/h. Here Ic_i is the indicator function which takes the value 1 in C-i and 0 
elsewhere. 

(c) If the event corresponds to propensity function a* and p{x,t) > 1/ft. for all x G C-i then 
place a particle in Ci. Remove a particle’s worth of mass from the PDE solution in the region 
C-i i.e. p{x, t + t) = p{x, t) — Ic_i/ft- 

(d) update the current time, t = tc- 
(viii) lftp<tc the the PDE regime is updated: 

(a) Update the PDE solution according to the numerical method described in Section 2.1 using 
p{x, t) as the previous value of the solution. 

(b) Update the current time, t = tp and set the time for the next PDE update step to be tp = 
tp -\- dtp. 

(ix) If t < tfinal, return to step (iii). 

Else end. 


Table 2. Algorithm for the coupling of the compartment-based regime with the PDE-based regime in 
the case of high copy numbers. 


to employ our algorithm. However, it is important to show that our algorithm is robust to such situations. 
In physical terms each of these simulations corresponds to particles, artificially kept in one half of a box, 
being allowed to diffuse into the other half of the box. 
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t=0 t=50 t=50 



XXX 


(a) (b) (c) 

Figure 4. Maintenance of (a) an initially uniform steady state by (b) Algorithm 1 and (c) Algorithm 
2. In each sub figure the blue curve is the solution of the PDE in flp, the blue histograms represent the 
concentration of particles in each of the compartments in Qq and the red curve represents the analytic 
solution of the mean-field model. Note that the blue curve is difficult to see since, for the most-part, it 
lies underneath the analytic solution. The number of particles is = 1000. The simulation results are 
averaged over 100 repeats. 


The mean-field model corresponding to this situation can be formulated as an initial-boundary-value 


problem: 





[-1,1], 

(19) 

with zero-flux boundary conditions 

dp 

dx 

= 0, 

(20) 

and initial conditions 





p{x, 0) = H{—x) for X G [—1,1], 

(21) 

or 





p{x, 0) = H{x) for X € [—1,1], 

(22) 


depending on which side of the box we wish to initialise the particles. This equation can be solved using 
separation of variables and Fourier series. For initial condition (21) the solution is given by: 


p{x,t) 


1 

2 


OO „ . 

^ (2n - 1)^ V 



(2n — 


(23) 


For initial condition (22) the minus in front of the sum is replaced with a plus. 

In Figure 5 we compare the analytically derived density to the density produced using our algorithms 
given initial condition (21). We see a favourable comparison between the solution according to our 
algorithms and the mean-field solution. A comparison of the densities through time for Algorithm 1 and 
2 can be seen in Movies M3 and M4 of the SM, respectively. In Movies M5 and M6 of the SM we present 
the same comparison but for a single realisation of Algorithms 1 and 2 respectively. In Figure 1 of the 
SM we compare the particle density of the two algorithms against the analytic solution of the mean-field 
model for initial condition (22). The comparison is again favourable. 
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Figure 5. Particles initially in the left-hand half of the domain (a) spread diffusively (with diffusion 
coefficient D = 1/400) throughout the domain eventually reaching steady state. By t = 1000 the system 
is assumed to have reached steady state. Panels (b) and (c) show the particle density profile according 
to Algorithm 1 at times 250 and 1000 respectively. Panels (d) and (e) show the same comparison for 
Algorithm 2. Figure descriptions and particle numbers are as in Fig 4. The simulation results are 
averaged over 100 repeats. 


3.3 Test Problem 3: A morphogen gradient formation model 


In this final test problem we model the formation of a morphogen gradient by considering basic first order 
reactions. Morphogen is produced at cc = — 1 with rate DX. We implement this production as a flux 
boundary condition in the PDF region of the model: 


dp 

dx 



= -A. 


(24) 


We also allow morphogen to decay uniformly across the domain with rate p and we implement a zero-flux 
boundary condition at a; = 1. These zeroth order (production) and first order (degradation) reactions 
allow for the description of the mean dynamics of the fully stochastic model by a PDF: 


dt 


= D 


d'^p 

dx"^ 


fip. 


(25) 


This is the mean-field PDF which we expect the density to satisfy across the whole domain. The analytic 
solution for the expected evolution of the density for this system (with boundary condition analogous to 
equation (24)) can therefore be found explicitly as 


p{x, t) = pst{x) + pu{x, t), 


( 26 ) 
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where 


Pst{,x) = A< 


ID 


cosh (^y^p,/D{l — x)'j 


^ cosh {2yJp,/D 

gives the steady state of the morphogen gradient and 

, , ^ cos(n7r(a; + l)/2) 

exp(-Ait) + 2^ ^ ^ ’ 


, . \D . 

Pu{x,t) = tanh I 2 a 


1 + D{mT)‘^/Ap 


exp - 


D{mr)‘^ 


+ p ] t 


(27) 


(28) 


Since equation (25) is the PDE we would expect the mean-field density of the fully stochastic model to 
satisfy, we use this as the PDE we solve in Dp in our hybrid model. In Dc we implement the usual 
diffusive ‘reactions’ and additional mono-molecular decay reactions in each compartment and a zero flux 
boundary condition at a; = 1. 

In Figure 6 we compare the analytically derived density to that produced using our algorithms given an 
initially empty domain. We increase the value of the diffusion coefficient by a factor of 10 (in comparison 
to previous simulations) to D = 10 x = 1/40 (in order to achieve a sensible steady state profile). 
We choose A = 1000 and p = 0.1, giving a rate of influx of particles into the domain of DX = 25. 
As in previous examples, we see a favourable comparison between the solution according to our hybrid 
algorithms and the expected mean-field solution given by equations (26)-(28). A comparison of the 
densities through time for Algorithms 1 and 2 can be seen in Movies M7 and M8 of the SM, respectively. 

A more quantitative study of the similarities between the solution of our algorithms and the analytic 
solutions is presented in the next section. 


4 Error analysis 


In this section we revisit the comparison between the mean particle densities predicted by our hybrid 
algorithm and the expected values given by the corresponding deterministic model. It is our aim in 
developing our algorithms that the error between the hybrid methods and the deterministic solution 
should be of the same order of magnitude as the accuracy associated with the individual numerical 
simulation techniques. In particular we developed our algorithms in order that the stochastic model 
should not ‘see’ (i.e. be influenced by) the interface. Since it is well established that the descriptions of 
particle density in Dp and Dc are valid representations of diffusion, we address the question of ‘whether 
our algorithms can correctly simulate the flux over the interface?’ In order to test this we compare 

NDc{t)= [ p{x,t)dx, (29) 

the expected number of particles in ^Ic the deterministic model to 

K 

NHc{t) = ^ (30) 


the expected number of particles in Dc in the hybrid models. Here Ai represents the mean number of 
particles in the compartment averaged over 100 repeats of the hybrid algorithms. For completeness 
we also compare 


Nopit) = J p{x,t)dx, (31) 

the expected number of particles in Dp in the deterministic model to 


NHp{t) = J p{x,t)dx, 


( 32 ) 



18 


(a) 




(b) (c) 



X X 


(d) 


(e) 


Figure 6. A morphogen gradient forms though a combination of influx at the left-hand boundary and 
decay of morphogen particles throughout the domain. Panel (a) shows the initially empty domain. By 
t = 60 the system is assumed to have reached an approximately steady state. Panels (b) and (c) show 
the particle density profile according to Algorithm 1 at times 15 and 60 respectively. Panels (d) and (e) 
show the same comparison for Algorithm 2. Figure descriptions are as in Fig 4. The simulation results 
are averaged over 100 repeats. 


the expected number of particles in flp in the hybrid models. 

The results for test problem 1 are shown in Figure 7. In simulations generated by our hybrid Algorithm 
2, the mass in the PDE-based regime is very similar to the expected mass given by the uniform solution 
of the diffusion equation. The same is true of the mass in the compartment-based regimes. Results are 
quantitatively indistinguishable (excepting for expected stochastic variation) for Algorithm 1, although 
for brevity we do not explicitly show the corresponding figures. 

In Figure 8 we display similar comparisons for test problem 2, in which the particles are initialised 
uniformly in flp only. As the simulation evolves we see that the expected mass predicted by the hybrid 
algorithm matches closely the expected mass given by the analytic solution in both ftp and flc- Further, 
the insets demonstrate that the relative error between the masses generated by the analytic solution and 
the hybrid model are low and have no systematic bias, but instead fluctuate about their expected values. 
Similar results hold for the initial condition in which the particles are initialised uniformly in 12(7 only 
(see Figure 2 of the SM). As in Figure 7, since the results of the simulations of our two hybrid algorithms 
are quantitatively similar, we present only the simulation results from our second hybrid algorithm. 

Finally in Figure 9 we display similar comparisons between the mass in each of the regimes for both of 
the hybrid algorithms we have developed when simulating the formation of a morphogen gradient (test 
problem 3). Again, a comparison of the results from the hybrid algorithm and the analytic solution shows 
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(a) (b) 


Figure 7. Error plots for the uniform particle distribution maintenance example. In panel (a) we 
display the normalised (i.e. divided by the total number of particles in the domain) expected value of 
the total mass in Qp given by the analytic solution (red, dashed line, Npip{t)/{Npip + Npcit))) and 
the normalised average mass in ilp according to our hybrid Algorithm 2 (black line, 

NHp{t)/{NHP + NHc{i)))- In panel 7(b) we display the same comparison but for the mass in the 
compartment-based regime (i.e. Np)c{i)/{NDP + Np,c{t)) (red, dashed line) and 
NHc{i)/{NHP + Npcit)) (black line)). The masses generated by our hybrid algorithm are averaged 
over 100 repeats as in previous density comparison figures. No systematic bias is discernible even for 
simulations run over longer time periods and with more particles. 


good agreement, especially for long times. However, there is a clear disparity between the mass in flc for 
Algorithm 2 and the expected mass in that region (see Figure 9 (d)). In particular the hybrid algorithm 
under-represents the expected mass in the compartment-based regime for early times. Although it is not 
clearly visible in Figure 9 (c) due to the scale of the axes, there is a corresponding over-representation 
of the expected mass in Up at early times for Algorithm 2. This is not the case for Algorithm 1 (see 
Figures 9 (a) and (b)). 

The reason for this disparity is a technical one. Recall that in Algorithm 2, in order to avoid negative 
solution values, a particle is not allowed to be removed from the pseudo-compartment, C-i, until the 
PDF solution across that compartment has value at least 1/h. This avoids the possibility that the removal 
of a pseudo particle may send the value of the PDF solution in C-i negative. As the domain is initially 
empty and fills from the PDF region there is a significant period of time in which mass should be flowing 
over the interface, but is restricted from doing so by an artificial barrier in our approximate Algorithm 
2. The flow of particles into Oc is, therefore, retarded in comparison to the analytic solution for which, 
of course, no such barrier exists. In Algorithm 1, however, we see no such problem (c/ Figure 9 (b)), 
because mass which is moved across the interface is removed from the whole of flp, rather than just the 
pseudo-compartment. This justifies our assertion that Algorithm 2 is suitable for simulations in which 
the number of particles at the interface are high, whereas Algorithm 1 is suitable for situations in which 
there may be low copy numbers at the interface. We explore the ramifications of these results further in 
the discussion. 

In order to determine the robustness of our coupling mechanism we performed further analysis of the 
error for test problems 1 and 2 (‘maintenance of a uniform steady’ state and ‘particles flowing from Up to 
flc’j respectively). Specifically we investigated the emergence of error associated with the discretisation 
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(a) (b) 

Figure 8. Error plots for test problem 2, in which all particles are initially in flp. In (a) we display the 
normalised expected value of the total mass in Up given by the analytic solution (red, dashed line, 
Ndp/{Ndp + Njjc) and the normalised average mass in flp (black line, Nhp/{Ndp + N^c) according 
to our second hybrid algorithm. In (b) we display the same comparison but for the compartment-based 
regime. In each panel, the inset displays the evolution of the relative error (i.e. {Nhp — Nj:,p)/Njjp or 
{Nhc — Np,c)/Np,c^ respectively) throughout the simulation. Again, the masses generated by our 
hybrid algorithm are averaged over 100 repeats. 


parameters dtp (the time discretisation associated with PDE updates) and h (the compartment size in 
flc) (see Figures 10 (a) and (b), respectively). The time-averaged relative error in the expected number of 
particles in He simulated using our hybrid simulation was determined by comparing to the same quantity 
predicted by a continuous mean-field analytical PDE solution to the respective test problems. In order 
to reduce the effect of stochastic fluctuation on the error measurement we averaged the particle numbers 
over 100 repeat simulations each each with N = 1000 particles for each of the test problems. Due to 
the computational intensity of these algorithms with such larger particle numbers and the wide variety 
of simulation parameters we conducted our error analysis using the algorithm two. We expect the error 
behaviour to be the same for algorithm 1 since the source of error we have identified (see below) should 
affect both algorithms in a similar manner. 

The results of varying both h and dtp separately are shown in Figure 10. On the vertical axis, error 
is defined as {{Npc — Nhc)/Ndc) where, as before, Np,c is the known expected number of particles 
in He and Npc is the averaged simulated number of particles in the same region. The angle brackets 
denote the time average of the relative error recorded at unit time intervals throughout the simulation. 
A positive error indicates a net bias in the simulations for particles towards Hp whereas a negative error 
indicates a net bias for particles towards He. 

Test problem 1 contains smaller sources of error than test problem 2. These errors have a consistent 
order of magnitude as that of expected stochastic fluctuations as a result of finite copy number. The 
difference in the magnitude of the error between test problem 1 and test problem 2 can be attributed to the 
way in which mass transfer is balanced or unbalanced. For test problem 1, net diffusive transfer of mass 
into the pseudo-compartment from the PDE is zero as is the net transfer of particles between the pseudo¬ 
compartment and He. In test problem 2 there is a net transfer of mass into the pseudo-compartment 
from the PDE regime and a net transfer out to the compartment-based regime. The systemic sources of 
error that we have identified in our method result from a non-zero net particle flux over the interface. 
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(a) (b) 



(c) 


(d) 


Figure 9. Error plots for the morphogen gradient formation example. Figure descriptions are as in 
Figure 8 


The pseudo-compartment method results in these errors being very small. 

Whilst still very small, a comparatively larger error associated with the hybrid modelling approach 
is observed in cases of non-zero net flux (in which the well-mixed assumption does not strictly hold 
within the pseudo-compartment). Increasing both h and dtp results in an increase in relative error in 
test problem 2 (see Figure 10). The reason for this increase is not directly related to the hybrid model 
but rather to the individual single-scale models used. Significant increases in the error for test problem 
2 can be seen when dtp is increased beyond 0{dx^/D). Similarly, an 0{h?) error is expected for large 
compartment sizes, h. In both cases, the effective diffusion is reduced (in comparison to that expected by 
the equivalent continuous PDF) as a result of discretisation error in the PDF and compartment region, 
respectively. The reduction of effective diffusion leads to a reduction in the net diffusive flux over the 
interface regardless of the implementation of the interface itself. The increase in error with increasing h 
and dtp is therefore, a manifestation of discretisation error in the simulation methods employed in flc 
and flp, respectively. 

The sign of the relative error in test problem 2 is negative for low h and low dtp. This is a manifestation 
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Figure 10. Error in compartment-based particle numbers varies with (a) PDE step size, dtp, and (b) 
compartment size, h, for both test problem 1 (blue) and test problem 2, in which all the particles are 
uniformly distributed in flp (red). Eor (a) simulation parameters are as given previously with dtp 
varying in steps of powers of ten from 10“"^ to 1 (note the log-scale on the horizontal axis.). For 
convenience in (b) we chose Ax = 1/256 and varied h in powers of two from 1/64 to 1/4. In test 
problem 1 final simulation time was 50 (as in Eigure 7), whereas in test problem 2 the final simulation 
time was 250. 


of a small increase in diffusive flux caused by our algorithm pseudo-compartment. The cause of this error 
is due to the breakdown of the assumption that mass in C-i is evenly distributed. Indeed, in test problem 
2, mass flowing into C-i from the remainder of Up is instantaneously considered to be evenly spread in 
C-i for the purposes of consideration of jumping into lie- Instead, this mass lies a small distance away 
from the expected position at the centre of the node on the side which is furthest from the interface. 
This erroneous evenly-distributed assumption results in a slightly increased net flux of particles over the 
interface that would otherwise be expected. The reason that this source of error is not prevalent in test 
problem 1 is because there is no net flux of particles flowing into or out of C-i that would contribute to 
such an error. The effect of this error is always to create a net increased bias over the interface in the 
direction of the net flow (an overall negative relative error in the case of test problem 2). It is important 
to notice that this error is common at the interface of hybrid models of this type. The inclusion of the 
pseudo-compartment rather than a clean interface forces particles which would contribute to this error 
to have a decreased number of transfer events associated with them and therefore compound a much 
smaller error (see for example, the two regime method Flegg et al. (2012)). 

Fine details of the error and its dependence on the small model parameters Ax, dtp and h is complex 
and difficult to determine due to other algorithmic factors (for example, particles need to wait between 
PDE times steps to be considered for C-i —>■ Ci when they are near the edge of C-i). It is also unclear 
if the reflecting boundary at I plays any role as a source of error for particular discretisation parameters. 
Some of these issues may be responsible for complex error behaviour in Figure 10 for small h and small 
dtp. On the other hand, details of fluctuation of the error for small h and small dtp are on the order of 
fluctuations that are expected as a result of finite copy number and may not be systemic in nature. 

Importantly, for small h and dtp, whilst there is a bias for particles to flow over the interface in the 
direction of net flow, this bias is still very small. The authors recommend the use parameters Aa:, dtp 
and h that individually minimize the discretisation error of the PDE and compartment-based methods 
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used. The interface should create an bias approximately proportional to h in the direction of net flow 
but this bias is significantly reduced due to the dual nature of pseudo-compartment. For completeness, 
the authors recommend that dxp << h to maintain the approximation of a PDF in flp although it would 
appear that the method is robust when this approximation is relaxed. Indeed, if dxp = h one would 
expect no interface-related error at all since transfer into the pseudo-compartment from the “PDF” is 
concentrated at the node at the centre of the pseudo-compartment. In this case though, the PDF becomes 
more closely described as a diffusion master equation on a lattice. 

5 Discussion 

We have introduced two algorithms which allow a coupling between a continuum PDF description of 
particle density in one region of the domain and a discrete stochastic compartment-based description in 
another. Using a pseudo-compartment in the PDF regime, which has properties of both descriptions, 
we were able to couple the two simulation methodologies together in such a way that the expected flux 
across the interface is maintained. We evidenced this though a series of examples, demonstrating a good 
correspondence between the expected mass in each of the modelling regions (averaged over several repeats 
of the hybrid simulation algorithms) and the mass in those regions determined by the corresponding fully- 
deterministic mean-field model. We also analysed the errors in the coupling method over a wide range 
of model parameter values and found them to be small (always less than 1%) even at the extremes of 
the parameter regimes in which the single-scale models we employ on either side of the interface are not 
expected to work well. 

Our method relies on a knowledge of the deterministic PDF description to which the mean of the 
compartment-based simulations corresponds. Whilst this is easy to obtain in the examples we have 
studied, when higher order reactions are introduced to the stochastic model there is no longer an exact 
closed-form representation of the expected particle density as a PDF in terms of the first moment alone. 
In order to obtain a self-contained equation we must make moment closure assumptions, which means that 
the PDF description will no longer correspond exactly to the mean behaviour of the compartment-based 
system. However, in situations where particle numbers are sufficiently large the closed PDF description 
will often be a close match to the expected densities of the compartment-based model and, as such, we 
expect that our algorithms will still be applicable. Furthermore, large particle numbers are precisely the 
regime in which we want to make use of the PDF description, further improving the applicability of our 
algorithms. 

It is precisely this consideration which will provide motivation for our next work on this subject to 
develop a hybrid method with an adaptive interface. In order to make the hybrid algorithm as efficient 
as possible, the interface should be able to move so that only regions with sufficiently low concentrations 
are represented by the compartment-based model and regions with high concentration are consistently 
represented with the PDF-based model. This would ultimately allow for Algorithm 2 to be used in all 
cases for efficiency reasons. A good set of rules that could be used to move an interface was explored in 
(Robinson et al., 2014). This will increase the computational efficiency of the algorithm, since we will 
no longer be simulating the individual particle movements in the discrete model when the continuum 
description can be tolerated. There is also scope to interface the algorithms we have developed in this 
manuscript with algorithms designed to couple compartment-based models and molecular-based models 
(Flegg et al., 2015, 2012) in order to build a “three-regime method”. Another important extension of this 
work is to consider compartment-based models on irregular lattices which will allow us to capture more 
complex and realistic domain geometries. 

For the sake of clarity in this manuscript, we presented the algorithms themselves and the illustrative 
examples in a single dimension, however, the algorithm is easily adapted for higher dimensional problems 
with (hyper)-planar interfaces. We therefore expect our algorithms to be highly applicable to stochastic 
simulations of systems in which molecular populations are high in some regions and low in others, meaning 
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that neither a fully-deterministic-continuum nor a fully-stochastic-individual-based model is appropriate. 
Our algorithms provide significant increases in speed in comparison to fully stochastic models and more 
accurate and representative modelling in comparison to fully-deterministic models. 
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1 Particles diffusing into an empty region from the compartment- 
based side 


Here we present the comparison between the analytical solution of the mean-field model and the density 
profile given by our algorithms for the initial-boundary-value problem from the main text: 


dt dx ^’ 


for X G [—1,1], 


( 1 ) 


with zero-flux boundary conditions 


dx 


= 0 , 

a;=0,l 


and initial condition 

p{x, 0) = H{x) for X G [—1,1]. 


( 2 ) 

( 3 ) 


This initial condition specifies that all the particles are initially in the compartment-based regime and 
none in the PDE-based regime. As such this simulation provides a good test of whether the algorithm 
is converting particles on the compartment-based side to mass on the PDE side appropriately. Clearly, 
because of the similarity between the mean-field solution and the density profiles given by our algorithms 
(see Figure 1) we can be confident that both of our algorithms are working correctly. In supplementary 
Movies M9 and MIO a comparison of the densities through time for algorithm 1 and 2, respectively, can 
be seen. 

In Figure 2 we compare the mass in each regime between our hybrid algorithm 2 and the analytical 
solution (given by equation (23) of the main text with the minus in front of the sum substituted for 
a plus) for the initial condition in which all particles are initialised in the compartment-based regime, 
He. The expected mass predicted by the hybrid algorithm matches closely the expected mass given by 
the analytical solution in both Hp and He. The insets demonstrate that the relative error between the 
masses generated by the analytical solution and the hybrid model are low and have no systematic bias, 
but instead fluctuate about their expected values in a similar manner to Figure 8 of the main text. 
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Figure 1. Particles initially in the right-hand half of the domain spread diffusively (with diffusion 
coefficient D = 1/400) throughout the domain eventually reaching steady state. Panel (a) gives the 
initial condition. By t = 1000 the system is assumed to have reached an approximately steady state. 
Panels (b) and (c) show the particle density profile according to algorithm 1 at times 250 and 1000 
respectively. Panels (d) and (e) show the same comparison for algorithm 2. Figure descriptions and 
particle numbers are as in Fig 4. The simulation results are averaged over 100 repeats. 
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(a) (b) 


Figure 2. Error plots for test problem 2, in which all particles are initially in ^c- In (n) we display the 
normalised expected value of the total mass in flp given by the analytical solution (red, dashed line) 
and the normalised average mass in flp (black line) according to our second hybrid algorithm. In 2(b) 
we display the same comparison but for the compartment-based regime. In each panel, the inset 
displays the evolution of the relative error (i.e. [Nhp — Ndp)/Ndp or {Nhc — Ndc)/Ndc^ 
respectively) throughout the simulation. Again, the masses generated by our hybrid algorithm are 
averaged over 100 repeats. 

















